Crystallization of classical multi-component plasmas 
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We develop a method for calculating the equilibrium properties of the liquid-solid phase transition 
in a classical, ideal, multi-component plasma. Our method is a semi-analytic calculation that relies 
on extending the accurate fitting formulae available for the one-, two-, and three-component plasmas 
to the case of a plasma with an arbitrary number of components. We compare our results to those 
of Horowitz, Berry, & Brown (Phys. Rev. E 75, 066101, 2007), who use a molecular dynamics 
simulation to study the chemical properties of a 17-species mixture relevant to the ocean-crust 
boundary of an accreting neutron star, at the point where half the mixture has solidified. Given 
the same initial composition as Horowitz et al., we are able to reproduce to good accuracy both 
the liquid and solid compositions at the half- freezing point; we find abundances for most species 
within 10% of the simulation values. Our method allows the phase diagram of complex mixtures 
to be explored more thoroughly than possible with numerical simulations. We briefly discuss the 
implications for the nature of the liquid-solid boundary in accreting neutron stars. 

PACS numbers: 89.90.+n, 97.60. Jd, 26.60.-c, 97.80.Jp 
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I. INTRODUCTION 



During the crystallization of a plasma containing mul- 
tiple ion species, the chemical composition of the solid is 
in general different from that of the liquid. This type of 
chemical separation is important for both white dwarfs 
[l[ and accreting neutron stars @- The interior of a 
white dwarf is a mixture of carbon, oxygen, and traces of 
other elements, most abundantly neon. As the star cools, 
chemical separation leads to the formation of an oxygen- 
and neon-rich core. The energy released through the 
gravitational settling of the denser core material heats 
the star and can delay cooling by several Gyr Q . A neu- 
tron star accretes mostly hydrogen and helium from its 
companion, but this material undergoes a series of nu- 
clear reactions, including rapid proton capture Q and 
then electron capture reactions y5|], to produce a variety 
of elements. Through accretion the mixture is pushed 
deep into the star and solidifies. Recent numerical simu- 
lations have shown that the mixture undergoes chemical 
separation during solidification [2|], possibly forming a 
two-phase solid [6j. The composition of the liquid ocean 
and the structure and composition of the crust have im- 
portant implications for a range of observed phenomena. 
For example, the resulting thermal conductivity deter- 
mines the cooling rate of transiently accreting neutron 
stars following extended accretion outbursts 7, 8]. The 
mechanical strength of the crust limits the size of a pos- 
sible crust quadrupole and therefore gravitational wave 
emission [9]. 

Several groups have studied the liquid-solid phase 
transition and chemical separation of two- and three- 
component plasmas in the classical, ideal limit (i.e., ig- 
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noring quantum mechanical effects on the ions and treat- 
ing the electrons as a uniform background; cf. Ref. fiof). 
Early works (e.g., Ref. ill]) studied phase transitions 
in carbon-oxygen plasmas, but the approximations used 
were too crude for application to the interiors of white 
dwarfs. Accurate calculations using the mean spherical 
approximation in the density-functional formalism were 
performed by Barrat et al. [Hj], who studied carbon- 
oxygen plasmas, and by Segretain and Chabrier who 
studied arbitrary two-component plasmas with atomic 
number Z ratios up to 2 (see also Ref. .14] , where carbon- 
oxygen- neon plasmas are examined). Using Monte Carlo 
calculations and Z ratios up to 5, Ogata et al. [H[ stud- 
ied arbitrary two- and three-component plasmas and Dc- 
Witt and Slattery [16| studied arbitrary two-component 
plasmas with a very accurate measurement of the liquid 
free energy (see also Refs. [IE EH)- All of these groups 
present phase diag rams as a function of ion abundance, 
and some [ill 1 1 6l ) also present fitting formulae for the 
liquid and solid free energies. Using these diagrams and 
fitting formulae, one can determine the phase transition 
properties for a two-component plasma of any ion type 
and abundance. 

These calculations are particularly useful for the inte- 
rior of a white dwarf, where there are only two or three 
dominant elements. But in the ocean of an accreting 
neutron star there are around 10-20 elements with abun- 
dances > 1% each one with a potentially important 
effect on the behavior of the phase transition and chem- 
ical separation of the mixture. The available analytic or 
numerical results for this type of system are extremely 
limited. We are aware of only one study of phase transi- 
tions in plasmas with more than three c omp onents, that 
of Horowitz et al. 0] (see also Refs. @, llSt). These au- 
thors used molecular dynamics simulations to study a 
17-component plasma with a composition similar to that 
expected at the ocean-crust interface of an accreting neu- 
tron star. Due to the large amount of computing power 
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necessary to run each simulation, the phase transition 
properties have so far only been calculated for one com- 
position. 

We present here a method for rapidly calculating the 
properties of the liquid-solid phase transition in a multi- 
component plasma in the classical ideal limit, for any 
initial composition and ion types. Our method is a semi- 
analytic calculation that relies on extending the accurate 
fitting formulae available for the one-, two-, and three- 
component plasmas to the case of a plasma with an arbi- 
trary number of components. We test our method using 
the one data point available for a plasma with more than 
three components, the calculation of Horowitz et al. Q , 
and show that it performs very well in that specific case. 

The paper is organized as follows. In Section |TT] 
we describe the semi-analytic calculation as it ap- 
plies to the one-component plasma (Section III Aj) . the 
two-component plasma (Section III B[) . and the multi- 
component plasma (Section III CI) . In Section IIIII we 
present our results for the 17-component mixture of 
Horowitz et al. Q. We conclude in Section [TV] The pres- 
sure term in the Gibbs free energy and its effect on the 
phase transition, the importance of the deviation from 
linear mixing for the liquid free energy, and a simplified 
derivation of the deviation from linear mixing for the 
solid free energy, are discussed in three appendices. 



II. METHOD 
A. The one-component plasma 

We assume in this paper that the system has reached 
equilibrium, i.e., the state of lowest free energy. The 
validity of this assumption and non-equilibrium effects 
such as diffusion and sedimentation will be discussed in 
a later paper. We also assume here that the phase transi- 
tion happens at constant volume, in which case the equi- 
librium configuration of the system is determined by the 
state with the lowest Hclmholtz free energy, F = U — TS. 
In reality the transition happens at constant pressure and 
at minimized Gibbs free energy. The error introduced by 
using the constant volume approximation is discussed in 
Appendix [SJ We find that for the mixture considered 
in Section IIIII the abundance in the liquid state of each 
ion species is in error by no more than 2%. While the 
percentage errors in the abundances in the solid state are 
typically larger by factors of ~ 2-5, the absolute errors 
for each ion species are similar in either state. (Since 
this trend holds true for most of the approximations we 
make in this paper, we hereafter quote errors in our ap- 
proximations only for the liquid abundances.) Note that 
in transitions at constant volume, the free energy of the 
electrons is identical in the liquid and in the solid and so 
has no effect on the properties of the phase transition. 

The Helmholtz free energy of the liquid or solid phase 
of a one-component plasma (OCP) can be described 
as a function of only the number of ions N, the tem- 



perature T, and the Coulomb coupling parameter T = 
{Ze) 2 /{ak B T) = Z 5 / 3 T e . Here Ze is the ion charge, a is 
the ion separation, and k B is the Boltzmann constant; 
T e = e 2 / (a e k B T) is the electron coupling parameter, 
where a e = [3/(47rn e )] 1 / 3 is the mean electron spacing 
and n e = ZN/V is the electron density. 

The ideal gas contribution to the free energy of a one- 
component plasma -Fideai is given by 
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where m t = Am p is the mass of the ion and (k B T) Ry = 
k B T2h 2 /(m l Z 4 e 4 ) is the thermal energy expressed in 
ionic Rydberg units. The free energy of the liquid 
phase of a one-component plasma -F/ OCP is well fit for 
Te [1,200] by 



pOCP 

/, OCP (D ee J, = - 0.899172r 



Nk B T 



1.8645r 32301 
0.2748 ln(r) - 1.4019. (2) 



The previous formula is from the Monte Carlo calcula- 
tions of DeWitt and Slattery [16[ , with the modification 
that the ideal gas contribution to the free energy [Eq. (JTJ] 
has been removed. Other formulae for fP CF can be found 
in Refs. QjJUJlIl] (see also Ref. [H,^); for the ran § e 
of r we are concerned with in this paper (15 < T < 200), 
the differences between these formulae, and between the 
numerical data these formulae are based on, are less than 
0.006. 

The free energy of the solid phase of a one-component 
plasma F° CP is well fit for V G [160, 2000] by 



^OCP 

Nk B T 



0.895929r+ 1.51n(r) - 1.1703 
10.84 176.4 5.980 x 10 4 
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r 3 



(3) 



The previous formula is from Dubin 25]; it was derived 
using a combination of analytic methods and a fit to the 
Monte Carlo calculations of Ref. [26j. As in the liquid 
case, we have modified Eq. ([3]) from its original form by 
removing the ideal gas contribution. Another formula 
for Fg F /(Nk B T) of similar accuracy (with less than 
0.004 difference from Ref. [25| or the numerical data for 
160 < r < 2000) can be obtained from the molecular 
dynamics calculations of Ref. [24| (see also Refs. [20l.[27T]). 
In this paper we neglect the T -2 and T~ 3 terms in Eq. © 
and use the following approximation for F® CP : 



fs(T) 



OCP 



Nk B T 



0.895929r+ 1.51n(r) 
10.84 



- 1.1703- 



(4) 
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This expression fits the numerical data for f60 < T < 
300 with an accuracy several times lower than that of 
Eq. © [differing by up to 0.02 for T ~ 160]. We use 
this expression in place of Eq. ©, however, because it 
behaves qualitatively better for small T, as we discuss 
below. 

The free energy difference based on these fits is given 

by 

Sfg CP (r)^ (/,-/ s ) OCP 

= - 0.003243r + 1.8645r ' 32301 

- 1.7748 ln(r)- 0.2316 + 10.84/T . (5) 

In equilibrium the system will be in the state of lowest 
free energy: when Sf OCP < 0, the OCP is in the liquid 
state, and when Sf OCP > 0, it is in the solid state. When 
fijOCP _ q t nere i s a phase transition between the liquid 
and solid state. This occurs at 



r crit = 178.6 



(0) 



in the above equation. Note that if we had used Eq. © 
instead of Eq. ([4]) to calculate <5/^ CP , we would obtain 
Tcrit = 175.3, which is in agreement with the most ac- 
curate estimates currently available for this value (e.g., 
r crit = 175.0 ± 0.4 in Ref. [Toj]); our r crit differs from the 
true transition value by about 2%. 

Equation © is only accurate for T e [160, 200]. While 
there are no Monte Carlo or molecular dynamics data 
available for /, OCP when T > 200, Ichimaru et al. [H| 
have calculated / ; OCP up to T — 1000 using the "im- 
proved hypernetted chain" (IHNC) method. For T 6 



[200,1000], if / ; UU1J is given by Ref. [28| and f ; 
given by Eq. ([4]), the approximation 

Sf OCP (p) = 0.09 + 0.0043(r - 200) 



OCP 



is 



(7) 



fits the free energy difference to within 0.2. This error is 
of similar magnitude to the error in the IHNC method 
for r > 200 (as extrapolated from comparisons between 
IHNC approximations and Monte Carlo calculations at 
T < 200; see, e.g., [ID, |3C|), and is several times smaller 
than the error that would be obtained by a direct appli- 
cation of Eq. © to the domain T G [200, 1000]. 

There are currently no published results (numerical 
or otherwise) for /; above T = 1000 or f® CP below 
r = 160. However, we expect 5f OCP to increase mono- 
tonically with increasing T, not just in [160, 1000] but for 
all T. In other words, for the OCP the solid state should 
always become more stable with respect to the liquid as T 
increases, and less stable as T decreases. Equation ([7]) ex- 
tended out to arbitrarily large T remains consistent with 
this assumption, but Eq. ([5]) extended down to T = 
does not. This is because Sf^ t CP decreases with T for 
T e [0,50]. An even stronger argument against (5/^ CP 
representing the true free energy difference at small P is 
that <5/fi t CP > for T < 17, which would imply that the 
OCP were in the solid state at very low T. Note that 
these effects are even worse if Eq. © is used to represent 



/° CP : in that case the free energy difference decreases 
with r for T 6 [0, 85] and is greater than zero for T < 51. 
To avoid small-I" problems, we cut off Eq. © at T = 100 
and assume that below this value the free energy differ- 
ence is given by 



(5/ UOP (r) = -0.37 + 0.0046(r - 100) 



(8) 



i.e., by the line tangent to Sf^ t CP at T = 100. If we had 
instead used Eq. © to represent / s OCP in 5fg t CP , Eq. © 
would change to <5/ OCP (r) = -0.30 + 0.0025(r - 100). 
Such a change leads to 'errors' in the multi-component 
results fSections lII Bl and lH C[) of no more than 5% for the 
liquid abundances, comparable to what is seen in Fig. [3] 
of Appendix [Bj 

Our final expression for Sf OCP , valid over all T, is 

p/° CP (r), ioo<r<20o, 
<s/ OCP (r) = < -0.37 + o.oo46(r - ioo) , r < ioo , 

I 0.09 + 0.0043(r - 200) , T> 200, 



where Sf° t CP (T) is given by Eq. ©. 



(9) 



B. The two-component plasma 

The free energy of a two-component plasma (TCP) can 
be described as a function of N, T, and the Coulomb cou- 

5/3 

pling parameter Tj = Z, t T e and fractional composition 
Xi = Ni/N of either species of ion. Here TV = TVi + 
is the total number of ions and n e = (Z\N\ + ^2^2)/^ 
is the total electron density. For the rest of this sec- 
tion we will identify the composition of the TCP by 
x\ and the Coulomb coupling parameter by Y\, since 
we can express xi and T2 as functions of these values: 
x 2 = 1 - xi and T 2 = {Z 2 /Z 1 ) 5/3 Ti. Note that through- 
out this paper we choose to label the ionic species such 
that Z\ < Z 2 < ■ ■ ■ < Z m , where m is the total number 
of species; Z\ always represents the ion with the smallest 
charge. 

The free energy of the liquid phase of a two-component 
plasma is given by 



/i TCP (F 1 ,x 1 )=E : 



i=i 



fP CP (r t ) 



In x 



^fl{Tx,xx), 



(10) 



where (Z) = Y)j—i XiZi is the average ion charge. The 
Y^i=i x i m i^i JzJ*) term is the (ideal gas) entropy of mix- 
ing for two species of volumes Z\Ni/n e and Z 2 N 2 /n e , 
and Afi is the deviation from linear mixing in the liquid. 
The deviation term A// has a similar dependence on Xi 
to the entropy of mixing term, but is in general much 
smaller in magnitude (see, e.g., Refs. [HI, [H, |3l| ) . We 
therefore expect this deviation to have a minimal effect 
on the phase transition properties for most systems. In 
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our calculation we set A/; = and use the linear mixing 
approximation: 



i=l 



(Z) 



(11) 



The error introduced by neglecting the A// term in the 
expression for /, TCP is discussed in Appendix [B] 

The free energy of the solid phase of a two-component 
plasma is given by 



/ s TCP (r l! ^)=E ; 



A/ S (r 1;2;i ), 



(12) 



where A/ s is the deviation from linear mixing in the solid. 
Unlike A/;, which is generally small even at large Ti 
(Appendix [5J , A/ s is comparable to the other terms in 
f s and grows linearly with we therefore expect Af s 
to play an important role in setting the phase transition 
properties. For charge ratios Rz — Z 2 /Zi in the range 
Rz €[1:5] the deviation is well fit by 



Af s {ri,Xi) ~T 1 x 1 x 2 Ag(x 2 ,Z 2 /Z 1 ) 



(13) 



where 

Ag{x,R z ) 



C(R 2 



1 + i+o'ife-'i) V^(V^ - 0.3)(Vi - 0.7)(Vi - 1) ' 

(14) 



C{R Z ) 



0.05(i?2 



[1 + 0.64(i? z - 1)][1 + 0.5{R Z - l) 2 ] 



(15) 



Equation (fTBj) is from the Monte Carlo calculations of 
Ogata et al. [l5|, and is accurate to within 10% for 
Rz 4.5; a similar formula (though accurate only for 
Rz < 2) can be found in DeWitt and Slattery To 
estimate the error introduced to our results by adopt- 
ing Eq. (|13|) . we run several calculations with a deviation 
of 1.1A/ S (ri,xi) and 0.9Af s (T u Xl ) [i.e., 10% higher or 
lower than the deviation we use in our model]. For the 
TCP, we find errors in the liquid abundances of 5% or 
less, with the largest errors at high T values and mod- 
erate charge ratios (Rz ~ 1.5). For the 17-component 
mixture and T value considered in Section HTT1 the errors 
in the liquid abundances are only 2% or less. 

For a TCP at a particular value of Ti, we find the 
state of lowest free energy as a function of composi- 
tion by u sing the "double-tangent" construction (see, 
e.g., Ref. [32(): We construct lines tangent to the min- 
imum free energy curve / m i n = min(/;,/ s ) in at least 
two points, corresponding to the compositions a\ and b\; 
an example of this construction is shown graphically in 
Fig. [U Any homogeneous composition x\ that lies be- 
tween eti and b±, i.e., any x\ which can be expressed as 




FIG. 1: (Color online) An example of the double-tangent con- 
struction, for R z = 34/8 and I\ = r crit /6 (cf. Figs. [4] and 
[5}. The stable compositions ai and 62 (i.e., 1 — ai and 1 — 61) 
are marked by filled circles; here, one of the mixtures is stable 
in the liquid state and one is stable in the solid state. Note 
that the curves /; and f s plotted in this figure are given not 
by Eq. (|lip and Eq. (|12[) . respectively, but by these equations 
minus the term Y^r=i x ifP CP (F*)- The values of 02 and 62 
obtained are the same whether /; TCP and fJ CI> or these mod- 
ified expressions are used: adding terms constant or linear in 
the Xi's to both free energy curves has no effect on the results 
of the double-tangent construction. 



Aa\ + (1 — A)b\ = X\ for some < A < 1, satisfies 
Af min (ai) + (1 - A)f min (bi) < fminixi) and is therefore 
unstable with respect to a heterogeneous mixture of a\ 
and b\. In this paper we refer to the locus of all points 
(Y\,x\) that lie between double-tangent points fFi,ai) 
and (Ti, 61) as the 'unstable region' of the phase diagram. 

Note that double-tangent points a\ and b\ can poten- 
tially be constructed from the liquid curve to itself, from 
the solid curve to itself, or from the liquid curve to the 
solid curve, depending on the behavior of /; and f s [see 
Eqs. (fTT|) and (fT2)) ]. In some cases 'triple-tangent' points 
can be constructed; typically this occurs when the solid 
curve is tangent to itself and to the liquid curve (when 
the liquid is at the "eutectic point"; see, e.g., Ref. [13|). 
The liquid-solid solutions are discussed below, in Sec- 
tion (|II B ip . In the approximation we have adopted 
above, where the deviation from linear mixing for the 
liquid is A/; = 0, tangents to the liquid curve fi do 
not intersect the curve at any other point [cf. Eq. (fTTj) ]; 
therefore there are no liquid-liquid solutions. Because of 
the A/ s > term in the solid curve, which grows pro- 
portional with Ti [see Eq. (fT3])]. when Ti is large enough 
there will always be regions of f s where double tangents 
can be constructed from the solid curve to itself. These 
solid-solid solutions will be examined in a later paper. 
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1. Solving for the liquid-solid equilibrium of the 
two- component plasma 

For a two-component plasma, liquid-solid phase tran- 
sitions occur at compositions and V values where double- 
tangent lines can be drawn connecting the free energy 
curves of the liquid and the solid. Under these condi- 
tions a liquid state of composition a\ and a solid state 
of composition 61 exist simultaneously as a mixture. For 
a double-tangent line connecting /; to f s , the line must 
satisfy 



//(a x ) = f' s (h) 



and 



(16) 
(17) 



fl(a 1 ) + (b 1 -a 1 )fl(a 1 )=f s (b 1 ). 
For later convenience we rewrite these equations as: 

/,(ox) + (1 - oi)//(o x ) = /.(61) + (1 - (18) 
and 

fifa) ~ 01/1(0!) = /.(61) - bifM ■ (19) 

Using Eqs. (JTTJ and (TT2|) the system of equations to solve 
becomes 



Zx 



r 1 b 2 l.Ag(b 2 ,Z 2 /Z 1 )-bxb 2 



dAg 



dx 



{b 2 ,Z 2 /Zx) 



(20) 



that to map out the full phase diagram we also need to 
know the shape of the solid-solid unstable region; this is 
most important at large Tx- Examples of phase diagrams 
for TCPs (including both types of unstable regions) are 
shown in the appendices. 



C. The multi-component plasma 

The free energy of an m-component plasma (MCP) can 
be described as a function of N, T, the fraction compo- 
sition of each ion species Xi — Ni/N (though x m is not 
needed, since x m = 1 — x i)i an d the Coulomb coupling 
parameter of one ion species. In the following we solve for 
Ti = zl ,S T e and then use the relation T; = {Z l /Z l f /3, T 1 
to find the other parameters. 

As with the two-component plasma, the free energy of 
the liquid phase of a multi-component plasma is very well 
described by the linear mixing rule (but see Appendix IBI): 



ff GV {Tx,xx,...,x m -x) 



zZ- 



fP CP (r t ) 



In 



(Z) 



(23) 



where (Z) = YaLx x i Z i- 

The free energy of the solid phase of the MCP is 



f^{Tx,x 1 ,...,x m . 1 ) 



zZ- 



Zi 



Af s (Fx,xx, ■ ■ -,x m -x) ■ 



(24) 



5/ OCP (r 2 ) + inU Z2 * " 2 



z 2 \ z 2 



V(Z)J (Z) b 
+ Txbx I Ag(b 2 ,Z 2 /Zi) - b x b 2 



dAg 



dx 



(b 2 ,Z 2 /Zx) 



(21) 



where (Z) a = ^£ J a i Z l , (Z)b = ^biZi, and 
dAg 



dx 



(x,Rz) = 

C{R z ){2x - 3^/5 + 1-21 - 0.105/Vi) 



1 + - 0-3)(Vi ~ 0.7)(Vi - 1) 

(22) 

[cf. Eq. (fT4l) ]. With these two equations (and ai+a 2 = 1, 
bi + b 2 = 1), if we are given Tx we can solve for ax and 
b\. This allows us to trace out the liquid-solid unsta- 
ble region of the phase diagram for Tx versus x\. Note 



According to Ogata et al. [15} . the deviation of the solid 
from linear mixing A/ s for a three-component plasma is 
given to good accuracy by 



A/ S (ri,xi, . . .,x m -x) 

rn — 1 m 

- ^ ^ TiXiXjAg 



x I x 7 Z ' ( 



Z 4 ) , (25) 



where Zx < Z 2 < ■ ■ ■ < Z m and Ag(x,Rz) is given by 
Eq. (TT4"|) . We assume here that Eq. (f2"5j) applies for all 
m > 2. A partial justification for this assumption is 
provided in Appendix [C] 

In the m-component plasma we construct (m — 1)- 
dimensional hyperplanes tangent to the minimum free 
energy surface in at least two points, corresponding to 
the compositions a and b. Any homogeneous composi- 
tion x that lies between a and b, i.e., any x which can 
be expressed as Aa + (1 — A)b — x for some < A < 1, 
is unstable with respect to a heterogeneous mixture of a 
and b. 



6 



1. Solving for the liquid-solid equilibrium of the 
multi- component plasma 

For a multi-component plasma, liquid-solid phase tran- 
sitions occur at compositions and T values where double- 
tangent hyperplanes can be drawn connecting the free 
energy surfaces of the liquid and the solid. For a double- 
tangent hyperplane connecting fi(a) to f s (b), the hyper- 
plane must satisfy 

$L(g) = $L(ft } »e[l,m-l] (26) 

Xjjjbl \AjJlJ 1 

and 

f l {a)+(b-a)-Vf l {a) = f s {b); (27) 



or 

//(«) + j|(a) -ci-Vm = f s (b) + ^-(b) - 6 - V/.(6) , 

i€[l,m-l] (28) 

and 

/, (a) - a • V// (a) - / s (6) - b ■ V/ s (6) (29) 



Using Eqs. ([23]) and (pM)) the system of equations to 
solve becomes 



-7—77- = In ( bi—^r 

(z) a V{z) b 



A/ S (ri,xi, . . . ,x m _i) 



v 6i + 6j Zj 



(6 i + 6 i ) 2 



dAg 

dx 

dAg 



h + bj'Zj 



dx 



bj + bj ' Zi 



(30) 



for i e [l,m]. Here (Z) Q = X a i^V (^)i> = Y,biZi, and 
[^] ( x > R z) is again given by Eq. ([H]). With these m 
equations (and y] a,j = 1, bj = 1), if we are given the 
liquid composition a we can solve for the solid compo- 
sition b and Coulomb parameter Ti at which the liquid 
and solid states are in equilibrium; if we are given b we 
can solve for a and IV In this manner we can trace out 
the liquid-solid unstable region of the phase diagram for 
Ti versus x. As in the TCP case, to map out the full 
phase diagram we also need to know the shape of the the 
solid-solid unstable region. 

Alternatively, if we are given an initial composition x 
and the fraction < A < 1 of the solution in the liquid 
state (or the fraction 1 — A in the solid state), we can 
solve for Ti and the compositions of both the liquid and 
solid mixtures in equilibrium. We have 2m — 1 unknowns, 
Oi,...,a m _i, 61, . . . , & m _i, and IV but in addition to 
the the m equations Eq. (|30l) above we have the m — 1 
equations 

Aa,i + (1 - A)bi = Xi , ie[l,m-l]. (31) 



III. RESULTS 

As described in Section HI Horowitz et al. @ [hereafter 
HBB] use a molecular dynamics simulation to study the 



phase transition of a 17-component plasma. A total of 
27, 648 ions are placed in a simulation volume of length 
727.5 fm on a side. At the start of the simulation 50% 
of the plasma is in the liquid state and 50% is in the 
solid state. There is a uniform composition throughout 
the volume, given by the results of Gupta et al. Q (who 
calculate the composition of an accreting neutron star at 
a density of 2 x 10 11 g/cm 3 , after the accreted material 
has undergone proton and electron capture and various 
other reactions). As the system evolves, the tempera- 
ture is adjusted so that approximately half of the plasma 
remains in the liquid state and half remains in the solid 
state. After a simulation time of 5x 10 6 fm/c, the simula- 
tion is run at constant energy until the total time reaches 
151 x 10 6 fm/c. The results of the numerical simulation 
are shown in Table HI The final temperature of the simu- 
lation is expressed in terms of T± as well as the 'average' 
Coulomb coupling parameter, T = (Z 5 / 3 )r e . For each 
entry in Table HI a statistical (y/Ni) error is provided. 

We have applied our semi-analytic calculation (Sec- 
tion []TC| to the same 17-component mixture as is con- 
sidered by HBB. In Eq. (|3T|) we set x to the 'initial' com- 
position given in Table HI and choose A = 0.5, such that 
we are solving for the equilibrium state where 50% of the 
mixture is liquid and 50% is solid. We then use Eqs. (f3"0")) 
and (|31l) to find the final composition of the liquid and 
solid states, a and b. The result is given in Table [TT1 For 
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TABLE I: Abundance of chemical element Z, for various mixtures from the numerical simulation of Horowitz et al. Q|. 
Abundances are provided for the initial mixture (in the column labeled 'Initial') and the final liquid and solid mixtures (in the 
columns labeled "Liquid" and "Solid", respectively). For each final mixture, the average charge (Z) and Coulomb coupling 
parameter F = {Z B ^ 3 )F e are provided as well. The percentage error for each entry is given by lOO/x/M, where Ni = XiN and 
N = 27, 648. 

HBB results 



(Z)i = 28.04, (Z) s = 30.48 
r a = 27.7, T ; = 233, T s = 261 



z 


Initial 


Liquid 


% Error 


Solid 


% Error 


8 


0.0301 


0.0529 


3 


0.0087 


6 


10 


0.0116 


0.0205 


4 


0.0021 


13 


12 


0023 

V 7 • \J\J 


0.0043 


9 


0006 


24 


14 


0.0023 


0.0043 


9 


0.0005 


27 


15 


0.0023 


0.0043 


9 


0.0004 


30 


20 


0.0046 


0.0055 


8 


0.0029 


11 


22 


0.0810 


0.1024 


2 


0.0616 


2 


24 


0.0718 


0.0816 


2 


0.0635 


2 


26 


0.1019 


0.1065 


2 


0.1017 


2 


27 


0.0023 


0.0025 


12 


0.0027 


12 


28 


0.0764 


0.0744 


2 


0.0746 


2 


30 


0.0856 


0.0773 


2 


0.0949 


20 


32 


0.0116 


0.0099 


6 


0.0130 


5 


33 


0.1250 


0.1079 


2 


0.1388 


1 


34 


0.3866 


0.3408 


1 


0.4297 


0.9 


36 


0.0023 


0.0012 


17 


0.0030 


11 


47 


0.0023 


0.0030 


11 


0.0013 


17 



each entry in Table llTl an error is provided in terms of the 
percent difference from the corresponding HBB result. 

The results of Table |TT] are relevant under equilibrium 
conditions, which in the accreting neutron star means 
that the particles solidify and diffuse through the liquid 
and the solid faster than new material is accreted. Here 
we attempt to estimate the importance of the diffusion 
rate on the overall results. In order to do that, we repeat 
our calculation done with 'instantaneous diffusion' (Ta- 
ble [TTJ) , this time assuming 'no diffusion' in the solid [36( ■ 
As in the equilibrium case, the calculation starts with the 
plasma in the liquid state with initial composition given 
by HBB, and ends when 50% of the plasma is liquid and 
50% is solid. Unlike in the equilibrium case, however, we 
solve Eqs. (|30]) and (j31~j) many times, each time produc- 
ing a small amount of solid material (1 — A -C 1). Solid 
particles created in one step are removed from considera- 
tion in all future steps, since we are assuming that these 
particles do not mix. The liquid composition (a) calcu- 
lated in one step is used as the 'initial' composition (x) 
in the next step. 

While an exact treatment of the 'no diffusion' limit 
would require solving Eqs. ([50)1 and (pn| on a particle-to- 
particle basis, we find that a good approximation can be 
obtained using 500 steps with Ak = 1 — 1/(1001 — k) for 
each step k. [The difference between the final abundances 
calculated using 50 steps with Ak = 1 — 1/(101 — k) and 
500 steps with Ak = 1 — 1/(1001 — fc), e.g., is less than 
0.2%.] The result is given in Table [TTTl Note that for 
this choice for Ak, the number of solid particles created 
is the same in each step. The average solid composition 



is given by 

50 

W = 7^E^' (32) 
k=i 

where b k is the composition of the solid particles created 
in the kth step. 

A comparison of Tables [H] and Mil shows that calcula- 
tions done under the two diffusion limits give very sim- 
ilar results. For example, the abundance differences be- 
tween these two calculations are generally much smaller 
than between either calculation and the results of HBB. 
Therefore, we conclude that the error introduced into our 
calculation by assuming instantaneous diffusion rather 
than the actual diffusion rate (whatever that may be) 
is small. Note that even though the rate of diffusion 
has very little effect on the average composition in the 
solid, it has a strong effect on the how that composition 
varies locally. For sufficiently low diffusion rates, lamel- 
lar sheets or other structures may form in the solid (see, 
e.g., Ref. [Hj]); these structures can have a strong effect 
on the thermal conductivity and strength of the crust. 

A comparison of Tables [TTI and Hill to Table Q] shows that 
the semi-analytic calculation does quite well at reproduc- 
ing the results of the HBB numerical simulation. All of 
the abundances from the semi-analytic calculation are 
with 65% of the HBB values, and most are significantly 
closer. Also, many of the table entries that match poorly 
between the two works correspond to chemical elements 
with very low abundances, i.e., those elements that are 
most affected by the finite size of the simulation. For 
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TABLE II: Abundance of chemical element Z, for the liquid and solid mixtures from our equilibrium calculation. Here, 
instantaneous diffusion is assumed (see text). For each mixture, the average charge (Z) and Coulomb coupling parameter 
r = (z 5 / 3 >r e are provided as well. The initial liquid mixture is given by its value from HBB, and the system is evolved until 
there is 50% liquid material, 50% solid material. The percent error for each entry is given by 100 x (entry — HBB)/HBB. 

Instant diffusion 
(Z)i = 27.667, (Z)s = 30.930 





Ti = 26.57, T ; 


= 218.3, r s = 


256.1 (Li error: -4%) 






z 


Initial 


Liquid 


% Error 


Solid 


% Error 


8 


0.0301 


0.0513 


-3 


0.0089 


+3 


10 


0.0116 


0.0197 


-4 


0.0035 


+66 


12 


0.0023 


0.0039 


-8 


0.0007 


+10 


14 


0.0023 


0.0040 


-7 


0.0006 


+22 


15 


0.0023 


0.0040 


-7 


0.0006 


+54 


20 


0.0046 


0.0073 


+32 


0.0019 


-33 


22 


0.0810 


0.1213 


+18 


0.0407 


-34 


24 


0.0718 


0.0947 


+16 


0.0489 


-23 


26 


0.1019 


0.1161 


+9 


0.0877 


-14 


27 


0.0023 


0.0024 


-4 


0.0022 


-19 


28 


0.0764 


0.0758 


+2 


0.0770 


+3 


30 


0.0856 


0.0759 


-2 


0.0953 


+0.5 


32 


0.0116 


0.0095 


-4 


0.0137 


+5 


33 


0.1250 


0.1013 


-6 


0.1487 


+7 


34 


0.3866 


0.3076 


-10 


0.4656 


+8 


36 


0.0023 


0.0018 


-12 


0.0028 


-8 


47 


0.0023 


0.0033 


+9 


0.0013 


+2 



TABLE III: As in Table ILTl except that diffusion is assumed to be negligible in the solid (see text). 

No diffusion 



(Z)i = 27.370, {Z) s = 30.680 
ri = 27.38, T ; = 221.2, Y a = 260.6 (ri error: -1%) 



z 


Initial 


Liquid 


% Error 


Solid 


% Error 


8 


0.0301 


0.0526 


-0.5 


0.0076 


-13 


10 


0.0116 


0.0204 


-0.5 


0.0028 


+34 


12 


0.0023 


0.0041 


-5 


0.0005 


-14 


14 


0.0023 


0.0041 


-4 


0.0005 


-5 


15 


0.0023 


0.0041 


-4 


0.0005 


+19 


20 


0.0046 


0.0077 


+39 


0.0015 


-47 


22 


0.0810 


0.1289 


+26 


0.0331 


-46 


24 


0.0718 


0.1018 


+25 


0.0418 


-34 


26 


0.1019 


0.1240 


+16 


0.0798 


-22 


27 


0.0023 


0.0025 


+1 


0.0021 


-23 


28 


0.0764 


0.0786 


+6 


0.0742 


-0.5 


30 


0.0856 


0.0753 


-3 


0.0959 


+1 


32 


0.0116 


0.0091 


-8 


0.0141 


+8 


33 


0.1250 


0.0953 


-12 


0.1547 


+11 


34 


0.3866 


0.2863 


-16 


0.4869 


+13 


36 


0.0023 


0.0017 


-18 


0.0029 


-4 


47 


0.0023 


0.0034 


+15 


0.0012 


-11 



example, the two entries that match the worst between 
Tables HI and ILTl the solid abundances of elements Z = 10 
and 15, are represented in the simulation by only 58 and 
11 ions, respectively. 

Figures [5] and [3J provide further comparison of our re- 
sults and those of HBB. Figure [2] (cf. Fig. 2 of HBB) 
presents in graphical form the data from Tables HI and HT1 
i.e., the final compositions of the liquid and solid states 
for both the HBB numerical simulation and our semi- 



analytic calculation. Figure [3J (cf. Fig. 6 of HBB) shows 
the ratio of the solid abundance to the liquid abundance 
versus atomic number Z for both works. 

Also plotted in Fig. [3J are the abundance ratios in the 
'two-component' approximation. In this approximation, 
the abundance ratios for each element are calculated as- 
suming the plasma is composed of only two ion species, 
the element itself and the most abundant element in the 
mixture (i.e., i — 15 or Z = 34; see Table UJ. The 
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initial composition of the mixture is chosen such that 
the ratio of the abundances of the two elements is the 
same as in HBB (e.g., xx/xis = 0.0301/0.3866, but now 
x\ + X\5 = 1); however, the results do not change much 
qualitatively if we choose some other scheme. As with 
the 17-component plasma, we solve for the point where 
half of the plasma is liquid and half is solid. Note that 
the Z = 34 abundance ratio is not plotted in Fig. [3] for 
this approximation, as its value is different for each two- 
element pairing. The two-component approximation re- 
produces the abundance ratio trend of the 17-component 
plasma, including the relatively constant behavior at low 
Z and the peak at Z = 34. It does not give accurate 
absolute values of the ratios, particularly for Z around 
Z = 34 (where the true solid-to-liquid ratio is greater 
than unity). 

The abundances listed in Table Hare the compositions 
of the HBB liquid and solid states at the end of the simu- 
lation. These results may not represent the true equilib- 
rium state of the mixture because of the finite run time of 
the simulation. To show this effect, the HBB abundance 
ratios are plotted in Fig. [3] using one of three symbols: 
for a given chemical element, if at the end of the simula- 
tion run the ratio is evolving upward in time, it is plotted 
with an upward-pointing triangle; if the ratio is evolving 
downward in time, it is plotted with a downward-pointing 
triangle; and if the ratio is not changing or is oscillat- 
ing upward and downward, it is plotted with a diamond. 
The determination of the evolution direction for each el- 
ement is made using data from the simulation time steps 
t 6 = i/(10 6 fm/c) = 71, 113, and 151, i.e., the last three 
time steps shown in Fig. 6 of HBB. If the abundance ra- 
tio decreases (increases) from = 71 to 113 and from 
t 6 = 113 to 151, and the total decrease (increase) across 
both time intervals is more than 0.1, the ratio is said to 
be evolving downward (upward) in time; otherwise the 
ratio is said to be stable. Note that, for the most part, 
the HBB results are evolving toward the equilibrium val- 
ues found in our calculation; this behavior is especially 
apparent for Z £ [20, 34], which is also where the abun- 
dance ratios differ in the two works by their largest values 
[37| . This suggests that the errors given in Tables [TT] and 
IIIII are strong upper limits to the actual accuracy of our 
calculation. 



IV. DISCUSSION 

Using results from simulations of one-, two-, and three- 
component plasmas, we have developed a method for 
calculating the equilibrium properties of the liquid-solid 
phase transition in a plasma with an arbitrary number 
of components, in the approximation of a classical ion 
plasma in a uniform electron background. We used this 
method to calculate the phase transition properties for a 
17-component plasma with a composition similar to that 
which might exist in the ocean of an accreting neutron 
star, and compared the results to those of a molecular dy- 
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FIG. 2: (Color online) Abundances 1 as a function of chem- 
ical element Z, for the final liquid and solid mixtures. Both 
the values from our equilibrium calculation ("Liquid" and 
"Solid", the large open squares and circles, respectively) and 
from the numerical simulation of HBB ("HBB", the small 
filled squares and circles) are shown. 
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FIG. 3: (Color online) Ratio of the solid abundance to the 
liquid abundance x s /xi as a function of chemical element Z. 
Both the values from our equilibrium calculation ("Current 
work", the open squares) and from the numerical simulation 
of HBB ( "HBB" , the filled diamonds and triangles) are shown, 
as are the values predicted from the 'two-component' approx- 
imation ("TCP approx.", the open circles; see text). If for 
a given element the HBB ratio is still evolving at the end 
of the simulation, it is plotted with a triangle that points in 
the direction of evolution; if the ratio is not changing or is 
oscillating up and down, it is plotted with a diamond. 



namics simulation done at the same composition (HBB 
0). We found that our method accurately reproduces 
the results of the HBB simulation. Two sources of error 
in the simulation may mean that our results represent 
the actual system even more accurately than this com- 
parison suggests: First, the finite size of the simulation 
introduces statistical errors which for some components 
are larger than the discrepancies between the two works. 
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Second, the system is still evolving at the end of the sim- 
ulation, with many components approaching the values 
predicted by our calculation. 

As in the simulation of HBB, we have followed the 
17-component mixture until it reaches the state of 50% 
liquid and 50% solid. Under these conditions, the term 
representing the deviation from the linear mixing rule for 
the solid, A/ s , is a perturbation on the other terms in 
the free energy of the solid [see Eq. (f24|)]. In principle our 
calculation can continue to larger fractions of solid, i.e., 
larger values of the Coulomb coupling parameter T. How- 
ever, because A/ s increases linearly with T and eventu- 
ally dominates the free energy, the calculation at T above 
the half-freezing point is more sensitive to the form cho- 
sen for Af s . There is some numerical confirmation of 
our simple approximation for A/ s , Eqs. ([T4")) and ([25]) . 
for two- and three-component mixtures at large I\ but 
only for a very limited set of parameters (see Ref. }15|). 
Further numerical simulations are necessary to test the 
validity of these equations at large T for general param- 
eters and (to > 3)-component plasmas. 

Another consequence of the large and positive A/ s 
term is that for certain compositions, it is energetically 
favorable for a single solid phase to separate into two or 
more solid phases (see Section IIIBj) . Such a phase sepa- 
ration occurs at large T in the 17-component plasma sim- 
ulated by Horowitz et al. [19j . With our calculation we 
have not yet found any two-solid mixtures that represent 
the lowest energy state of the HBB plasma, in part be- 
cause the shape of free energy surface for the solid phase 
is very complicated at large T. We leave a more careful 
study of the solid-solid unstable region for future work. 

Once these issues are resolved, our calculation will al- 
low the complete phase diagram of multi-component mix- 
tures to be determined. We expect that these results 
will have important implications for the structure of the 
liquid-solid boundary in accreting neutron stars. For ex- 
ample, for an ocean temperature of T = 10 8 Ts K, an 
O-Se mixture with the same proportion of oxygen and 
selenium as in the HBB mixture (i.e., ~ 10%-90%) will 
begin to freeze at a density of p ~ 2x 10 7 Tg (/i e /2) g/cm 3 , 
where p e is the mean molecular weight per electron. As- 
suming that accretion is slow enough that the liquid and 
solid can come into equilibrium at each depth, our phase 
diagram for a charge ratio Rz — 34/8 in Fig. @] (or Fig. [5]) 
shows that the mixture will reach 50% solid within a fac- 
tor of two in density, but that complete freezing will not 
occur until much deeper, by a factor of ~ (34/8) 5 ~ 1400 
in density (corresponding to p ~ 3 x 10 10 Tg g/cm 3 ). 
This is a very different picture than the sharp transition 
between liquid and solid expected for a one-component 
plasma, and assumed in previous work on accreting neu- 
tron stars. Further work is needed to understand the 
effects of the various time-dependent processes that are 
active concurrent with accretion in the ocean-crust tran- 
sition layer, such as crystallization, diffusion, and sed- 
imentation. For example, sedimentation of the heavier 
solid particles could be important at low accretion rates, 



narrowing the transition layer. 
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Appendix A: The Helmholtz free energy versus the 
Gibbs free energy 

Because phase transitions in stars occur at constant 
pressure, not constant volume, the energy which is at a 
minimum when the system is in equilibrium is the Gibbs 
free energy, i.e., G = F + PV . We discuss here how our 
results (Section IIII[) change when the Gibbs free energy, 
rather than the Helmholtz free energy, is used to deter- 
mine the equilibrium state. 

To calculate the Gibbs free energy, we follow the per- 
turbation method of Ogata et al. [15| , though we ignore 
terms due to the electron exchange energy (see, e.g., 
Ref. [33j]; these terms are small for highly- relativistic 
plasmas such as are found at the ocean-crust bound- 
aries of accreting neutron stars). In the degenerate in- 
teriors of white dwarfs and neutron stars, the electrons 
make the dominant contribution to the total pressure 
(Pi ~ a(Z 2 / 3 )P e for r > 1; see, e.g., Ref. [H), and so 
we can treat the ion partial pressures as perturbations. 

The Helmholtz free energy of the system is 



F — Fg + F\ 



(Al) 



where Fg is the kinetic energy of the electrons and F\ is 
the free energy of the ions (the electron exchange term is 
ignored and the Coulomb term is folded into the ion free 
energy). The total pressure of the system is 



dF 



dF Q dFt 



p - —- w (A2) 



dV 



Let Vo be the volume of the unperturbed system, when 
only electrons contribute to the total pressure; let V01 be 
the volume of the perturbed system, when both ions and 
electrons contribute to the total pressure. Then the total 
pressure can also be expressed as 



and 



P = -Fo(V ) 



SV 2 



P = P(V ) + SVP'(V Q ) + —P"(V ) 



(A3) 



(A4) 



where SV — Vbi — Vq and we are using the notation 

p'(v ) = [&] v=Vo , 



etc. From Eqs. (|A2)| - (|A4|) . and 
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assuming SV is small (which can easily be checked a pos- 
teriori), we obtain 



SV 



F{(V ) 



The Gibbs free energy can be written as 

G = G(V ) + 5VG'(V ) + ^G"(V a ) + ■ ■ 

= F q (Vq) + F 1 (V a ) + P(V )V Q + SVP'(V )V 

SV 2 SV 2 
+ d ^-P"(V )V + d -^P'(Vo) + --- 



Fo(V ) + F^Vo) + PV 



[F{(Vo)f 



(A5) 

(A6) 

(A7) 
(A8) 



where in going from Eq. (|A6|) to Eq. (|A7|) we have made 
use of the thermodynamic relation 



V = - 



dG 
dP 



(A9) 



The Gibbs free energy is obtained from Eq. (|A8I) , once 
the value of Vq is known. For a given total pressure P, 
the volume Vq is determined by Eq. (|A3I) : We have (e.g., 
Ref. M) 



P 



8tt 2 A^ 



1 



In 



where the "relativity parameter 
A C (3^V) 1/3 - 



y 



PF 

m e c 



9tt\ 1/3 k B T 



am e r 



- y 

(A10) 



(All) 



is evaluated at V — Vq. Here a = e 2 /(Tic) is the fine 
structure constant and A c = h/(m e c) is the reduced 
Compton wavelength. The volume Vq depends only on 
the total pressure of the system, and so is the same for 
both the liquid and solid states. The Helmholtz free en- 
ergy in the unperturbed state, Fo(Vq), is also the same 
for both states. We can therefore ignore the F (Vq) an d 
PVq terms in Eq. (|A8[) when calculating the state of low- 
est free energy. Using 



F o(V ) 



1 m P C 2 



y 



V Q 9n 2 Xi y/T+y 2 



(A12) 



we arrive at our final expression for the Gibbs free energy 
of the liquid (i — I) or solid (i — s) state: 



9i 



G, 



Nk B T 
/i(T e ) - 



3(18tt) 1 /3(Z) y 



dfi 



(A13) 



where fi is the Helmholtz free energy given in Sec- 
tions EaUrn y( P) is found from Eq. (|AT0|) . and r e (y) 
is found from Eq. (|All|) (i.e., T e is evaluated at V = Vq)- 
We calculate the phase diagrams for two-component 
plasmas with charge ratios Rz = Z%/Z\ up to 34/8, first 
using the relevant expressions for fi and f s from Sec- 
tion [TIB] (i.e., ignoring pressure terms), and then using 
Eq. (|A13I) (including pressure terms). Note that the T e 
values in Eq. (|A13[) are evaluated at V = Vq , while those 
in Section III Bl are evaluated at ~ Vqi . In order to show 
the two sets of phase diagrams on the same axis we use 
the relation [cf. Eq. (|A5j)]: 



r e (y i) =T e (V ) 



r e (v ) 



Vqi 

Vq 



1/3 



2a dfi 



1/3 



(A14) 



where all instances of y and T e on the right-hand side 
of Eq. (|A14I) are evaluated at V — Vq. Here we choose 
to solve for r e (Vbi) of the liquid, although the results 
are practically the same if r e (Voi) of the solid is used 
instead (since the two r e values differ by at most 0.004% 
even for Rz — 4). Our results, plotted as a function of 

ri(Vbi) = Z 1 5/3 r e (V r i), are shown in Fig. H Not sur- 
prisingly, we obtain results very similar to those found 
by [la]: the assumption of transitions at constant vol- 
ume rather than at constant pressure has no effect on 
the phase diagram unless Rz ^ 2, in which case the 
only effect is to widen the unstable region slightly. For 
2 < Rz < 5 the unstable region widens by at most 1-2%, 
with the largest change occurring for Ti < r cr ;t. Since 
the calculation of Section Hill was done at a relatively low 
value of r (at Tz=% — 27, which is below r cr ; t for all 
species Z < 25), we expect that the results shown there 
will not change when the Gibbs free energy is used. At 
large T, however, when nearly all of the mixture is in the 
solid state (see Section IIV|) , inclusion of the Gibbs free 
energy in the equations of Section III CI may be necessary 
to accurately determine the phase transition properties 
under these conditions. 



Appendix B: The deviation from linear mixing in 
the liquid 

In our calculation we assume perfect linear mixing in 
the liquid state, by setting A// = 0. We discuss here how 
our results (Section Mil) change when a more accurate 
form for A// is used. 

There are several fitting formulae of A/; available in 
the literature (e.g, Refs. [lE Ell HI1 ) ■ We choose to im- 
plement the fit from Equation (9) of Potekhin et al. [35[ 
(hereafter PCCDR), since it provides accurate results for 
Afi over a wide range of T values, Z ratios, and fractional 
abundances of each species. It is also the only fit we are 
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FIG. 4: (Color online) Phase diagrams for charge ratios Rz = 
34/20 (top panel) and Rz = 34/8 (bottom panel). Phase 
transitions at constant volume are labeled "Helmholtz" , and 
transitions at constant pressure are labeled "Gibbs". To 
maintain consistency with earlier works (e.g., Refs. [l3l. [l5|). 
r^ 1 in units of T~ r ] t is plotted versus X2, where Z2 = 34 for 
all transitions. The unstable regions are marked by dots. The 
mixture is liquid for (x2, Pcrit/ri) points entirely above the 
unstable region; for points below any part of the unstable re- 
gion (such as the peninsula in the bottom-left corner of the 
top panel and the banana-shaped island in the bottom panel) 
the mixture is solid. 



aware of that is immediately applicable to plasmas with 
more than two components, though we do not make use 
of that feature here. 

We calculate the phase diagrams for two-component 
plasmas with charge ratios Rz up to 34/8, first for Afi = 
0, and then using Eq. (9) of PCCDR (i.e., for Afi ^ 
0). Our results are shown in Fig. [5J We find that the 
assumption Afi = has no effect on the phase diagram 
unless Rz It 3, in which case the only effect is to shift 
the I0W-X2 side (the left side, in Fig. [5j of the unstable 
region toward even smaller values of x^.- The shift is 
most significant for large Rz and T, with shifts of around 
5% of the width of the unstable region for Rz — 4 and 



FIG. 5: (Color online) Phase diagrams for charge ratios 
Rz = 34/26 (top panel) and Rz = 34/8 (bottom panel). 
Phase transitions where the liquid deviation term Afi is ig- 
nored are labeled "A/; = 0" , and transitions where the liquid 
deviation is given by Eq. (9) of Potekhin et al. [35| are labeled 
"PCCDR". 



Ti ~ r cr i t . Since our calculation was done at a relatively 
low value of T, we expect that the results of Section HTT1 
will not change when an accurate form for Afi is used 
(cf. Section . At larger values of T, a Afi term may 
be necessary to ensure the accuracy of the calculation. 

Here and in Section [X] we have compared phase di- 
agrams generated by our calculation to those that are 
generated if additional terms are considered. We can 
also compare our phase diagrams to those of other works. 
Particular fruitful comparisons can be made with Segre- 
tain and Chabrier [l3| and Ogata et al. [lEj . since these 
works present phase diagrams at several different values 
of Rz] the Rz values in Figs. 2] and [5J were chosen in 
part because of the similarity to the ratios presented 
in these two works (i.e., R z = 34/26 ~ 4/3 = 1/0.75, 
Rz = 34/20 ~ 5/3 ~ 1/0.55, and R z = 34/8 ~ 13/3). 
Our diagrams agree closely with those of [15[, with one 
important exception: for most values of Rz, this group 
finds 'azeotropic points' or eutectic points at X2 ^ 0.04 
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that do not exist in our diagrams. The close agreement 
for x 2 > 0.04 is due to the fact that both our group 
and theirs used fitting formulae with the same form for 
Af s [Eq. ([13] , while the poor agreement at x 2 < 0.04 is 
due to the fact that we used A// = while [15[ used a 
form for A/; that was negative for x 2 ^ 0.05. Our di- 
agrams agree less closely with those of [lj], though the 
agreement is still very good at small T (in the upper 
half of each diagram). Even at large T the diagrams of 
our group and theirs are qualitatively similar, with the 
main differences being the larger amount of stable solid 
regions at high x 2 and the delayed (in terms of increas- 
ing Rz) transition from spindle type to azeotropic type 
in the diagrams of [l3| . We find that the transition from 
spindle- type to azeotropic-type phase diagrams occurs at 
Rz — 1.2 ~ 28/34 ~ 1/0.83, which is a somewhat lower 
value of Rz than found by Segretain and Chabrier [l]| 
or DeWitt et al. [3 (1/0.72 ~ 1.4). 

Appendix C: The deviation from linear mixing in 
the solid 

In this section we provide a simple estimate of A/ s 
for multi-component plasmas, using the approximation 
that only nearest neighbors contribute to the interaction 
energy of each ion (see, e.g., Ref. |32j). The expression 
found here is too simplistic for use in our calculation, 
but illustrates the general form of A/ s for plasmas with 
three or more components; the A/ s term of Section |TT] 
[Eq. (US])] has a very similar form. 

Let Uij = Uij I (NksT) be the interaction energy be- 
tween nearest-neighbor ions of species i and j (uij — Uji). 
When all ion species are completely separated, the inter- 
action energy per ion for species i is uu/2, and the total 
interaction energy of the system is given by 

U S cp = ~Z ^ XiUa . (CI) 



When the ion species are mixed, the interaction energy 
per ion for species i is J2j x j u ij/%> assuming that the var- 
ious ions are randomly distributed throughout the mix- 
ture. The total energy of the system is then 

^mix ^ ^ ^ ^ ^ ^i^jUij . (02) 
i 3 

The internal energy of mixing for the solid, Au s — u m ; x — 
itscp, is given by 



i V i I 

i j>i ' 

The free energy of mixing can be found from the ther- 
modynamic identity 




where /3 = l/(fcsT) (see, e.g., Ref. Assuming that 

the interaction energies «y scale linearly with /? (which 
is true, e.g., if Uij oc r e ), we have 

\r- XX **** ("v - , (C5) 

i j>i 

which of the same form as Eq. ([231) . 
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